home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Dr. Windows 3
/
dr win3.zip
/
dr win3
/
PROGRAMR
/
GSRC208A.ZIP
/
MATRIX.C
< prev
next >
Wrap
C/C++ Source or Header
|
1993-08-26
|
6KB
|
272 lines
#include "copyleft.h"
/*
GEPASI - a simulator of metabolic pathways and other dynamical systems
Copyright (C) 1989, 1992 Pedro Mendes
*/
/*************************************/
/* */
/* matrix operations */
/* */
/* MICROSOFT C 6.00 */
/* QuickC/WIN 1.0 */
/* ULTRIX cc */
/* GNU gcc */
/* */
/* (include here compilers that */
/* compiled GEPASI successfully) */
/* */
/*************************************/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include "globals.h"
#include "globvar.h"
#define NO_MTX_ERROR 0
#define NO_MULT 1
#define NON_INVERT 2
#define OUT_OF_MEM 3
#define NEARLY_ZERO 1e-10
#define TINY 1.0e-20;
/* reset a vector to zero */
void zerovct( double m1[MAX_MET] )
{
register int i;
for(i=0; i<MAX_MET; i++)
m1[i] = 0.0;
}
/* add two vectors */
void addvct( double m1[MAX_MET], double m2[MAX_MET],
double m3[MAX_MET], int r )
{
register int i;
for(i=0; i<r; i++)
m3[i] = m1[i] + m2[i];
}
/* add two matrixes */
void addmtx( double m1[], double m2[], double m3[], int r, int c )
{
register int i,d;
d = r * c;
for(i=0; i<d; i++) m3[i]=m1[i]+m2[i];
}
/* add one matrix with a scalar (double) multiple of another */
void addmtxsm( double m1[MAX_MET][MAX_MET], double m2[MAX_MET][MAX_MET],
double s, double m3[MAX_MET][MAX_MET],
int r, int c )
{
register int i,j;
for(i=0; i<r; i++)
for(j=0;j<c;j++)
m3[i][j] = m1[i][j] + s * m2[i][j];
}
/* add one vector with a scalar (double) multiple of another */
void addvctsm( double m1[MAX_MET], double m2[MAX_MET],
double s, double m3[MAX_MET],
int r )
{
register int i;
for(i=0; i<r; i++)
m3[i] = m1[i] + s * m2[i];
}
/* (s1 * v1) + (s2 * v2) - vectors */
void addvctsm2( double m1[MAX_MET], double s1,
double m2[MAX_MET], double s2,
double m3[MAX_MET], int r )
{
register int i;
for(i=0; i<r; i++)
m3[i] = s1 * m1[i] + s2 * m2[i];
}
/* multiply one matrix by another */
void multmtx( double m1[MAX_STEP][MAX_MET], float m2[MAX_MET][MAX_MET],
double m3[MAX_STEP][MAX_MET], int r1, int c1, int r2, int c2 )
{
register int j,k;
int i;
double acum;
/* if (c1!=r2) return(NO_MULT); disabled for speed */
for(i=0;i<r1;i++)
for(j=0;j<c2;j++)
{
acum=0.0;
for(k=0;k<c1;k++) acum += m1[i][k] * (double) m2[k][j];
m3[i][j]=acum;
}
/* return(NO_MTX_ERROR); disabled for speed */
}
/* scalar (double) multiple of multiplication of one matrix by a vector */
void multmtxvsm( double m1[MAX_MET][MAX_MET], double m2[MAX_MET],
double s, double m3[MAX_MET],
int r1, int c1, int r2 )
{
register int i,k;
double acum;
/* if (c1!=r2) return(NO_MULT); disabled for speed */
for(i=0;i<r1;i++)
{
acum=0.0;
for(k=0;k<c1;k++) acum+=m1[i][k] * m2[k];
m3[i] = s * acum;
}
/* return(NO_MTX_ERROR); disabled for speed */
}
/* multiplication of one matrix by a vector */
void multmtxv( double m1[MAX_MET][MAX_MET], double m2[MAX_MET],
double m3[MAX_MET], int r1, int c1, int r2 )
{
register int i,k;
/* if (c1!=r2) return(NO_MULT); disabled for speed */
for(i=0;i<r1;i++)
{
m3[i] = 0.0;
for(k=0;k<c1;k++) m3[i] += m1[i][k] * m2[k];
}
/* return(NO_MTX_ERROR); disabled for speed */
}
/* mutiply a matrix by a scalar (double) */
void multmtxs( double m1[MAX_MET][MAX_MET], double s,
double m2[MAX_MET][MAX_MET], int r, int c )
{
register int i,j;
for(i=0;i<r;i++)
for(j=0;j<c;j++)
m2[i][j] = m1[i][j] * s;
}
void back_sub(double (*a)[MAX_MET][MAX_MET], int n, int *indx, double *b)
{
int i, ii, ip, j;
double sum;
for( i=0, ii=-1; i<n; i++ )
{
ip = indx[i];
sum = b[ip];
b[ip] = b[i];
if( ii != -1 )
for( j=ii; j<i; j++ ) sum -= (*a)[i][j] * b[j];
else if( sum ) ii = i;
b[i] = sum;
}
for( i=n-1; i>=0; i-- )
{
sum = b[i];
for( j=i+1; j<n; j++ )
sum -= (*a)[i][j] * b[j];
b[i] = sum/(*a)[i][i];
}
}
int lu_decomp(double (*a)[MAX_MET][MAX_MET], int n, int *indx, double *d)
{
int i, imax, j, k;
double big, dum, sum, temp;
double vv[MAX_MET];
*d = (double) 1;
for( i=0; i<n; i++ )
{
big = (double) 0;
for( j=0; j<n; j++ )
if( ( temp = fabs( (*a)[i][j] ) ) > big ) big = temp;
if( big == (double) 0 )
{
return NON_INVERT;
}
vv[i] = (double) 1/big;
}
for( j=0; j<n; j++ )
{
for( i=0; i<j; i++ )
{
sum = (*a)[i][j];
for( k=0; k<i; k++ ) sum -= (*a)[i][k] * (*a)[k][j];
(*a)[i][j] = sum;
}
big = (double) 0;
for( i=j; i<n; i++ )
{
sum = (*a)[i][j];
for( k=0; k<j; k++ )
sum -= (*a)[i][k] * (*a)[k][j];
(*a)[i][j] = sum;
if( ( dum = vv[i]*fabs( sum ) ) >= big )
{
big = dum;
imax = i;
}
}
if( j != imax )
{
for( k=0; k<n; k++ )
{
dum = (*a)[imax][k];
(*a)[imax][k] = (*a)[j][k];
(*a)[j][k] = dum;
}
*d = -(*d);
vv[imax] = vv[j];
}
indx[j] = imax;
if( (*a)[j][j] == (double) 0 ) (*a)[j][j] = TINY;
if( j != n-1 )
{
dum = (double) 1/((*a)[j][j]);
for( i=j+1; i<n; i++ ) (*a)[i][j] *= dum;
}
}
}
int lu_inverse( double (*m1)[MAX_MET][MAX_MET], double (*m2)[MAX_MET][MAX_MET], int n )
{
double dt;
register int i,j;
double tvv[MAX_MET];
int emptyrow, mat_indx[MAX_MET]; /* indexes for matrix row permut. */
for( i=0; i<n; i++)
{
emptyrow = 1;
for( j=0; j<n; j++)
if ( (*m1)[i][j] != (double) 0 ) emptyrow = 0;
if ( emptyrow ) return NON_INVERT; /* if one row is empty m1 singular */
}
if ( !( j = lu_decomp(m1, n, mat_indx, &dt) ) ) return j;
for( j=0; j<n; j++ )
{
for( i=0; i<n; i++ )
tvv[i] = i==j ? (double) 1 : (double) 0;
back_sub(m1, n, mat_indx, tvv);
for( i=0; i<n; i++ ) (*m2)[i][j] = tvv[i];
}
return NO_MTX_ERROR;
}